In this model with $k$ variables and $d$ categories, let $x = (1,x_1, x_2,\ldots, x_k)$ be one data point including a constant $1$ for the intercept. There are $d-1$ column vectors $\beta_2, \beta_3, \ldots, \beta_d,$ each of length $k+1,$ that give the relative chances of each category in the ratios

$$p_1:p_2:\cdots:p_d = 1:e^{x\beta_2}:e^{x\beta_3}:\cdots:e^{x\beta_d}.$$

The random response for this data point places a specified number $s$ balls into $d$ bins (one for each category) according to these probabilities.

Thus, to simulate such data you need to

Specify the variables--their number $k,$ the number of data points $n,$ and all their values.

Specify the number of categories and the $\beta_j$ parameters.

Specify the value of $s$ for each data point.

Compute the probabilities determined by (1) and (2) according to this model.

Use those sizes (3) and probabilities (4) to generate random multinomial outcomes.

This leads to a straightforward R implementation. It randomly generates and stores the variables (1) in an array X and, given a randomly-generated array of coefficients (2) and randomly-generated size array (3), computes the probabilities (4) and applies rmultinom to each data point (5) to obtain the matrix of responses (one column per category) and store it in the array y.